function [tides_out] = compare_obsROMS_tides(basename,romsgaugedir,obsdir,gaugeNames,constit)
% compare_obsROMS_tides.m
%
% [tides_out] = compare_obsROMS_tides(basename,romsgaugedir,obsdir,gaugeNames,constit)
%
% Performs a tidal analysis of one year's record of sea-level from both
% observations and model output at multiple stations and compares the
% results statistically.  The tidal analysis is done using the T_Tide
% package, but is only run if both observations and a model extraction
% exist for each site.  Stations can be a combination of US and Canada
% sites, and the tidal constituents to compare can be specified.  Based on
% previous code from S. Giddings that was based on code from D. Sutherland.
%
% A 'switch' option exists to make simple comparison plots of the low-pass
% filtered sea-level time series.  To generate better plots of the model-
% observational comparisons, use 'plot_compare_obsROMS_tides.m' instead
% (which calls this function).
%
% 
% USES:
%   Processed hourly sea-level data from both the US and Canada stored in
%   .mat files
%
%     CAN_****_TIDE_yyyy.mat      Canada data (**** = station ID, yyyy = year)
%     NOAA_*******_TIDE_yyyy.mat  US data     (******* = station ID, yyyy = year)
%
%   Hourly tide gauge extractions from the numerical model 
%     'basename'_'station'.mat    (e.g., T40final_base_2005_noramp_LaPush.mat)  
%
%
% REQUIRES:
%     T_Tide tidal analysis package (http://www.eos.ubc.ca/~rich/index.html)
%
%
% INPUT:  (optional, can be specifed in code below)
%
%     'basename'      -  model run basename
%     'romsgaugedir'  -  path location of ROMS tide gauge extractions
%     'obsdir'        -  path location of observational tide gauge records
%                        currently assumes the following structure:
%                          obsdir/NOAA/NOAA_year
%                          obsdir/Canada/CAN_year
%     'gaugeNames'    -  vector of strings of tide gauge station names
%                        (e.g., 'La Push','Seattle,'Vancouver', etc.)
%     'constit'       -  vector of tidal constituent names for which to
%                        perform the analysis, e.g., {'M2';'S2';'K1';'O1'}
%
%
% OUTPUT:
%     'tides_out'     -  data structure with lots of info including:
%
%      tides_out.name          -  names of the tidal constituents
%      tides_out.freq          -  frequencies for the tidal constituents
%      tides_out.station_name  -  names of tide gauge stations analyzed
%      tides_out.lat           -  latitudes of tide gauges
%      tides_out.lon           -  longitudes of tide gauges
%      tides_out.WS            -  Willmott skill scores for the hourly
%                                 model-observation sea-level comparisons
%                                 at each tide gauge site
%      tides_out.WSlp          -  Willmott skill scores for the low-pass
%                                 filtered model-observation sea-level
%                                 comparisons at each tide gauge site
%      tides_out.SS            -  skill score statistics for the hourly
%                                 model-observation sea-level comparisons
%                                 at each tide gauge site
%      tides_out.SSlp          -  skill score statistics for the low-pass
%                                 filtered model-observation sea-level
%                                 comparisons at each tide gauge site
%
%     both SS and SSlp have fields such as:
%        tides_out.SS.SS       -  skill scores
%        tides_out.SS.CC2      -  correlation coefficients (r^2)
%        tides_out.SS.var      -  variance
%        tides_out.SS.mean     -  means
%     (refer to skill_score.m for more information)
%
%     and data structures for each constituent such as:
%      tides_out.M2
%
%     having fields
%        tides_out.M2.amp   -  amplitudes for the constituent at each site
%        tides_out.M2.pha   -  phases for the constituent at each site
%        tides_out.M2.ampE  -  amplitude errors for the constituent at each site
%        tides_out.M2.phaE  -  phase errors for the constituent at each site
%
%
% R. McCabe 05/2014

     
%% Set global paths and directories
addpath('../../alpha/');
Tdir = toolstart;
addpath(Tdir.tools);
addpath(Tdir.data);
addpath(genpath([Tdir.tools 'validation_tools']));
cwd = pwd;


%% User changeable zone... switches, flags, etc.
% Switches
plotsamplefigure = 0; % 1 = make simple plots of low-pass sea-level at each station

% Set default path locations, basename, gaugeNames, constit if none were provided
%
% Define which tidal constituents to compute
if ~exist('constit','var')
    constit = {'M2';'S2';'K1';'O1'}; %;'n2';'q1'};
end
%
% Set the run basename
if ~exist('basename','var')
    basename = 'T40final_base_2005_noramp';
end
%
% Location of ROMS tide gauge extractions
if ~exist('romsgaugedir','var')
    romsgaugedir = '/data1/rmccabe/runs/PostProcess/ROMS_tides/';
end
%
% Parent directory for the observations
if ~exist('obsdir','var')
    %obsdir = '/data1/rmccabe/tools_data/tide_gauge';
    obsdir = [Tdir.data 'obs_data/tide_gauge'];
end
%
% Set list of tide gauge locations to extract - comment out the ones you don't want!
if ~exist('gaugeNames','var')
    gaugeNames = {'Charleston'
        'South Beach'
        'Garibaldi'
        'Astoria'
        'Toke Point'
        %'Westport' % exists for 2006 onward
        'La Push'
        'Bamfield'
        'Tofino'
        'Neah Bay'
        'Port Angeles'
        'Victoria'
        'Port Townsend'
        'Friday Harbor'
        'Patricia Bay'
        'Cherry Point'
        'Vancouver'
        'Point Atkinson'
        %'Campbell River' % just outside of Cascadia domain
        'Seattle'
        'Tacoma'};
end


%% Location for the barometric pressure data files
% This is a temporary approach until a more accurate atmospheric correction
% is incorporated into the processed tide gauge data files themselves
metdir = [Tdir.data 'obs_data/ndbc/sg_baro/'];


%% Convert gaugeNames to a "no space" version for file names
gaugeNames_nospace = cell(size(gaugeNames));
for ii = 1:length(gaugeNames)
    gn = gaugeNames{ii};
    gn(gn == ' ') = [];
    gaugeNames_nospace{ii} = gn;
    clear gn
end
clear ii


%% Initialize tidal parameter storage for amplitudes and phases
for ii = 1:length(constit)
    constname = char(constit(ii)); % convert from cell to char
    % When we have model output to compare with then the first column is
    % for tide gauge observations, and the second column is for the model
    % tides_out.(constname).amp = NaN*ones(length(TG),2); obs = (:,1), model = (:,2)
    tides_out.(constname).amp = NaN.*ones(length(gaugeNames),2);
    tides_out.(constname).pha = NaN.*ones(length(gaugeNames),2);
    tides_out.(constname).ampE = NaN.*ones(length(gaugeNames),2); % amplitude error
    tides_out.(constname).phaE = NaN.*ones(length(gaugeNames),2); % phase error
end
clear ii constname


%% Loop through each ROMS tide gauge extraction...
disp(' ');
disp(' ***** Tide Gauge comparison ***** ')
disp(' ');
for ii = 1:length(gaugeNames)
    disp([' Station ',gaugeNames{ii},':  station ' num2str(ii) ' of ' num2str(length(gaugeNames))]);
    
    % If the ROMS file doesn't exist continue to the next loop iteration
    if ~exist([romsgaugedir basename '_' gaugeNames_nospace{ii} '.mat'],'file')
        disp([' Station ' gaugeNames{ii} ' was not found.  Moving on...'])
        continue;
    end
    
    % Load the extracted ROMS tide gauge file
    disp([' Loading ROMS extraction file ',basename,'_',gaugeNames_nospace{ii},'.mat']);
    load([romsgaugedir basename,'_',gaugeNames_nospace{ii},'.mat']);

    % Determine current year based on first ROMS output file
    year = datestr(ssh_moor.td(1),'yyyy');
    yyyy = str2double(year);
    td0 = datenum(str2double(datestr(ssh_moor.td(1),'yyyy')),1,1);
    
    
    %% Move to load in the corresponding observations
    % Get the unique station ID...
    [tgloc_id] = tidegauge_name2id(gaugeNames{ii});

    % Open the various NOAA folders looking for the matching obs .mat file...
    % Look for a NOAA file first
    cd([obsdir,'/NOAA/NOAA_',year]);
    fname = ['NOAA_',tgloc_id,'_TIDE_',year,'.mat'];
    if ~exist(fname,'file')
        % If the NOAA file didn't exist, look for a Canada file
        disp(['  File ',fname,' does not exist.  Trying Canada...']);
        cd([obsdir,'/Canada/CAN_',year]);
        fname = ['CAN_',tgloc_id,'_TIDE_',year,'.mat'];
        if ~exist(fname,'file')
            % If both the NOAA and Canada files no not exist, move on...
            disp(['  File ',fname,' does not exist.  Moving to next site...']);
            clear fname tgloc_id
            % Also clear out ROMS information
            clear ssh_moor roms_td roms_zeta
            continue
        end
    end

    % Load the OBS .mat file
    disp(['  Loading file ',fname,'...']);
    load(fname);
    clear tgloc_id fname
    
    % Move back
    cd(cwd);
    
    
    %% Tidal analysis of OBS and ROMS tide gauge records
    %
    % T_Tide analysis of observed sea-level
    [t_struct, ] = t_tide(tgauge.sea_level,'interval',1,'start time',...
        tgauge.td(1),'output','none','latitude',tgauge.lat,'error','wboot');
    
    cname = cellstr(t_struct.name); % cell array of all t_tide names
    for jj = 1:length(constit)
        constname = char(constit(jj));
        id = find(strcmpi(constname,cname)); % locate that particular constituent in t_tide output
        tides_out.(constname).amp(ii,1) = t_struct.tidecon(id,1);
        tides_out.(constname).ampE(ii,1) = t_struct.tidecon(id,2);
        tides_out.(constname).pha(ii,1) = t_struct.tidecon(id,3);
        tides_out.(constname).phaE(ii,1) = t_struct.tidecon(id,4);
        if (ii==1)
            tides_out.freq(jj) = t_struct.freq(id);
            tides_out.name(jj) = cname(id);
        end
    end
    clear jj cname constname id t_struct
    disp(' ');
    
    % Also store the station_name, lon, lat
    tides_out.station_name{ii} = tgauge.station_name;
    tides_out.lat(ii) = tgauge.lat;
    tides_out.lon(ii) = tgauge.lon;
    
    
    % T_Tide analysis of model sea-level
    [t_struct, ] = t_tide(ssh_moor.zeta,'interval',1,'start time',...
        ssh_moor.td(1),'output','none','latitude',ssh_moor.lat_rho,'error','wboot');
    %
    cname = cellstr(t_struct.name); % cell array of all t_tide names
    for jj = 1:length(constit)
        constname = char(constit(jj));
        id = find(strcmpi(constname,cname)); % locate that particular constituent in t_tide output
        tides_out.(constname).amp(ii,2) = t_struct.tidecon(id,1);
        tides_out.(constname).ampE(ii,2) = t_struct.tidecon(id,2);
        tides_out.(constname).pha(ii,2) = t_struct.tidecon(id,3);
        tides_out.(constname).phaE(ii,2) = t_struct.tidecon(id,4);
    end
    clear jj cname constname id t_struct
    disp(' ');
    
    
    %% Rename variables for convenience
    % Rename ROMS time and surface heights
    roms_td = ssh_moor.td; % should be a single row or column vector
    roms_td = roms_td(:)'; % make it a row vector
    %
    roms_zeta = ssh_moor.zeta;
    roms_zeta = roms_zeta(:)'; % row vector
    
    % Rename OBS time and surface heights
    obs_td = tgauge.td;
    obs_td = obs_td(:)'; % row vector
    %
    obs_zeta = tgauge.sea_level;
    obs_zeta = obs_zeta(:)'; % row vector
    
    
    %% Interpolate the observed gauge data to the ROMS times
    % This is for the skill calculations
    % NOTE:  both records should be hourly on the hour, but this makes sure
    % ALSO:  ROMS output generally starts at 01-Jan 01:00 instead of
    % midnight and extends to 01-Jan 00:00 of the following year instead of
    % 31-Dec 23:00...
    obs_zetai = interp1(obs_td,obs_zeta,roms_td);
    obs_tdi = roms_td;
    
    
    %% Inverse barometer correction
    % This is the approach adopted from S. Giddings and is approximate
    % ON THE TO-DO LIST FOR REVISION / IMPROVEMENT    
    switch gaugeNames{ii}
        case {'Charleston','South Beach','Garibaldi','Astoria'}
            NDBCname = 'nwpo3'; % Newport, OR, site
        otherwise
            NDBCname = 'ttiw1'; % Tatoosh Island, WA, site
    end
    [datem,~,pcorrectRaw] = InverseBarometerCorrection(metdir,yyyy,NDBCname);
    pcRaw = interp1(datem,pcorrectRaw,roms_td);
    obs_zetai = obs_zetai+pcRaw;

    
    %% Low-pass filter observations and model sea-level
    % Godin low-pass filter
    % NOTE: Need to include case for dealing with interior NaNs better
    obs_zetailp = Z_godin(obs_zetai');
    roms_zetalp = Z_godin(roms_zeta');
    
    % Flip back...
    obs_zetailp = obs_zetailp';
    roms_zetalp = roms_zetalp';
    
    
    %% Skill score and Willmott Skill score calculations
    % Recall input order is:  SS = skill_score(model, obs, etc); and WS = willmott_score(model, data);
    %
    % Skill calculations for de-meaned tidally-varying surface height
    tides_out.SS(ii) = skill_score(roms_zeta-nanmean(roms_zeta),obs_zetai-nanmean(obs_zetai),1);
    tides_out.WS(ii) = willmott_score(roms_zeta-nanmean(roms_zeta),obs_zetai-nanmean(obs_zetai));
    %
    % Skill calculations for de-meaned low-pass filtered surface height
    tides_out.SSlp(ii) = skill_score(roms_zetalp-nanmean(roms_zetalp),obs_zetailp-nanmean(obs_zetailp),1);
    tides_out.WSlp(ii) = willmott_score(roms_zetalp-nanmean(roms_zetalp),obs_zetailp-nanmean(obs_zetailp));
    
    
    %% Simple plots for low-pass filtered sea-level
    if plotsamplefigure
        
        figure; hold on;
        plot(roms_td-td0,obs_zetailp-nanmean(obs_zetailp),'b'); % observations
        plot(roms_td-td0,roms_zetalp-nanmean(roms_zetalp),'r'); % ROMS
        plot([0 365],[0 0],'k')
        legend('tide gauge','ROMS')
        xlim([0 365])
        
        title(['Low-pass filtered sea-level ' strrep(basename,'_','\_') ' at '...
            gaugeNames{ii} ', SS = ' num2str(tides_out.SSlp(ii).SS,2) ...
            ' r^2 = ' num2str(tides_out.SSlp(ii).CC2,2)]);

        %print(gcf,'-dpng',[outdir basename 'lp_' gaugeNames{gi}])
    end

    clear obs_td obs_tdi obs_zeta obs_zetai obs_zetailp
    clear roms_td roms_zeta roms_zetalp
    clear tgauge ssh_moor
end
clear ii count
disp(' ');


%% Print statistics of comparisons to the screen
disp([' ']);
disp(' ***** Tide Gauge comparison results ***** ')
disp([' ']);
disp([' Model run:  ',basename]);
disp(' Model-Obs sea-level comparison statistics:');
disp(' (SS = Skill score, r = correlation, WS = Willmott score)');
disp(' ');
for gi = 1:length(tides_out.lat)
    if gi == 1
        fprintf('\t \t hourly\t hourly\t hourly\t lp\t lp\t lp\n');
        fprintf('Station\t \t SS\t r\t WS\t SS\t r\t WS\n')
    end
    % To align screen output we need an extra tab if station names are less
    % than 8 characters long...
    if length(char(tides_out.station_name(gi))) < 8
        fprintf('%s\t\t %1.2f\t %1.2f\t %1.2f\t %1.2f\t %1.2f\t %1.2f\n',...
            char(tides_out.station_name(gi)),...
            tides_out.SS(gi).SS,sqrt(tides_out.SS(gi).CC2),tides_out.WS(gi),...
            tides_out.SSlp(gi).SS,sqrt(tides_out.SSlp(gi).CC2),tides_out.WSlp(gi));
    else
        fprintf('%s\t %1.2f\t %1.2f\t %1.2f\t %1.2f\t %1.2f\t %1.2f\n',...
            char(tides_out.station_name(gi)),...
            tides_out.SS(gi).SS,sqrt(tides_out.SS(gi).CC2),tides_out.WS(gi),...
            tides_out.SSlp(gi).SS,sqrt(tides_out.SSlp(gi).CC2),tides_out.WSlp(gi));
    end
end
disp([' ']);

